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Abstract 

A  multiscale  modeling  approach  for  the  prediction  of  material  stiffness  of  the  ionic  polymer  Nafion 
is  presented.  Traditional  rotational  isomeric  state  theory  is  applied  in  combination  with  a  Monte 
Carlo  methodology  to  develop  a  simulation  model  of  the  conformation  of  Nafion  polymer  chains 
on  a  nanoscopic  level  from  which  a  large  number  of  end-to-end  chain  lengths  are  generated.  The 
probability  density  function  of  end-to-end  distances  is  then  estimated  and  used  as  an  input  pa¬ 
rameter  to  enhance  existing  energetics-based  macroscale  models  of  ionic  polymer  behavior.  Several 
methods  for  estimating  the  probability  density  function  are  compared,  including  estimation  using 
Johnson  distributions,  Bezier  distributions,  and  cubic  splines. 


1.  Introduction 


Ionic  polymers  have  received  significant  attention  in  the  last  decade  due  to  the  property  that  the 
electromechanical  coupling  produced  through  the  transport  of  charged  and  uncharged  ions  within 
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the  polymer  matrix  provides  them  with  unique  transduction  capabilities  for  a  variety  of  appli¬ 
cations.  Specifically,  mechanical  deformations  produce  an  internal  redistribution  of  charges  and 
solvent  which  produces  a  measurable  electrical  signal,  thus  giving  the  polymer  sensing  capabilities. 
Reciprocally,  applying  an  electric  field  to  the  structure  produces  ionic  migration  toward  the  negative 
surface  which,  when  combined  with  the  redistribution  of  solvent,  leads  to  bending  in  the  structure, 
enabling  the  ionic  polymer  to  act  as  an  actuator.  As  an  actuator,  the  material  has  greater  actuation 
displacement  and  can  be  operated  under  much  lower  voltages  than  other  electroactive  materials 
(such  as  peizoelectric  ceramics).  These  materials  are  superior  to  a  number  of  existing  transducers 
for  several  reasons,  including  their  ability  to  operate  under  cryogenic  conditions  and  their  high 
gravimetric  energy  density.  The  growing  number  of  potential  applications  for  ionic  polymers  in¬ 
clude  artificial  robotic  muscles  and  biomimetic  sensors.  A  detailed  review  of  the  early  history  of 
these  materials  and  their  applications  is  provided  by  Shahinpoor  et  al.  [17]. 

One  factor  inhibiting  the  widespread  practical  implementation  of  ionic  polymers  is  the  fact  that 
many  aspects  of  their  behavior  are  not  fully  understood.  In  recent  years,  significant  focus  has  been 
directed  towards  the  development  of  suitable  models  which  quantify  their  unique  behavior.  Several 
models  have  focused  on  characterizing  the  electromechanical  response  through  the  first  principles 
of  fluid  transfer  [1,  4,  6,  16].  Other  models  represent  the  mechanical  deformation  of  ionic  polymers 
as  functions  of  current  [1,  13]  or  ion  transport  [11].  Most  recently,  there  has  been  progress  on 
formulation  of  models  which  incorporate  both  the  sensing  and  actuation  properties  of  the  material 
[14,  15], 

It  is  commonly  understood  that  the  polymer  chains  making  up  the  material  are  composed  of 
a  hydrophobic  backbone  with  side  chains  that  terminate  in  hydrophilic  ionic  groups.  In  the  bulk 
material,  these  polymer  chains  interact  with  each  other  so  that  the  bulk  material  has  at  least  two 
distinct  phases  —  a  hydrophobic  region  comprised  of  the  backbones  of  the  chains  and  a  hydrophilic 
region  comprised  of  the  ionic  terminal  groups.  When  the  material  is  solvated  with  water,  Hsu  and 
Geirke  [9]  propose  that  the  hydrophilic  ionic  terminal  groups  and  the  water  that  has  been  taken  up 
by  the  material  cluster  together.  The  model  further  suggests  an  idealized  structuring  whereby  the 
embedded  clusters  are  of  essentially  constant  radius,  uniformly  distributed  throughout  the  material, 
and  interconnected  by  channels. 
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In  Weiland  and  Leo  [23],  an  energetics  approach  is  used  to  study  the  equilibrium  state  of 
the  hydrophilic  cluster  regions  of  the  ionic  polymer  Nafion.  The  primary  goal  of  that  work  is  to 
better  understand  how  the  equilibrium  state  may  facilitate  or  inhibit  ion  transport  as  it  relates 
to  transduction.  One  significant  drawback  to  that  approach  is  the  lack  of  appropriate  nanoscale 
stiffness  data.  Moreover,  multiscale  stiffness  predictions  of  ionomers  in  general,  to  date,  have  not 
been  studied  and  represent  an  area  critical  to  the  tailoring  of  material  response. 

The  overall  goal  of  this  investigation  is  to  enhance  the  existing  energetics-based  models  of  ionic 
polymer  behavior  by  estimating  the  material  properties  of  the  hydrophobic  region  of  the  material. 
Specifically,  we  focus  on  the  development  of  a  multiscale  modeling  approach  for  the  prediction  of 
material  stiffness.  This  method  involves  the  application  of  rotational  isomeric  state  theory  (RIS), 
traditionally  used  for  studying  deformation  trends  and  material  properties  as  described  by  Flory  [8] , 
in  combination  with  a  Monte  Carlo  methodology  to  develop  a  simulation  model  which  characterizes 
the  conformation  of  Nafion  polymer  chains.  The  material  stiffness  depends  on  multiple  parameters 
including  the  effective  length  of  the  polymer  chains  comprising  the  material.  This  model  will  be  used 
to  simulate  a  large  number  of  distances  r  between  cluster  interaction  points  —  that  is,  points  where 
the  ionic  terminal  groups  for  one  pendant  chain  come  within  a  certain  distance  of  an  embedded 
cluster.  It  is  therefore  assumed  that  communication  between  an  ionic  terminal  group  and  a  cluster 
acts  as  a  crosslink.  For  the  single  Nafion  chain  in  Figure  1,  there  are  five  cluster  interaction  points, 
resulting  in  the  r  values  {r*  :  i  =  1, . . . ,  5}.  The  simulated  r- values  can  then  be  used  to  estimate 
a  probability  density  function  P[r )  for  cluster  interaction  distance  (or  equivalently,  polymer  chain 
length).  Our  modeling  methodology  is  based  on  the  approach  of  Mark  and  Curro  [12]  for  short 
chain  polymers  displaying  rubber  elasticity  and  is  similar  to  the  application  of  the  Mark-Curro 
method  to  study  the  effect  of  particle  reinforcement  in  polymers,  as  described  in  [18,  19]. 

Relation  of  Density  to  Macroscopic  Properties 

The  estimated  density  function  P(r )  for  polymer  chain  length  may  then  be  used  in  a  macroscopic- 
level  model  to  compute  Young’s  modulus  to  quantify  material  stiffness.  Whereas  details  of  the  de¬ 
velopment  of  the  corresponding  macroscopic-level  model  are  provided  in  [22],  in  brief,  Boltzmann’s 
approach  to  statistical  thermodynamics  is  used  to  relate  the  material  entropy  to  the  density  P(r) 
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Backbone 


Figure  1:  Distances  r  between  cluster  interaction  points  for  a  single  Nation  polymer  chain. 


through  the  expression 


S(r)  =  kin  P(r) 


where  k  is  Boltzmann’s  constant.  Under  the  assumption  of  rubberlike  elasticity,  the  “three  chain” 
model  as  described  by  Treloar  [20]  may  be  applied  to  compute  the  change  in  entropy  from  an 
unperturbed  configuration  to  a  distorted  configuration, 


AS 


S(rQa)  +  2S(r0a  l^2)  -  3S(r0) 


where  u  is  the  number  density  of  network  chains,  rQ  is  the  root  mean  square  of  the  simulation¬ 
generated  r  values  and  a  =  L/Li  is  the  relative  length  of  the  sample.  RIS  theory  assumes  that 
under  load,  the  rotation  about  a  given  bond  is  unrestricted  and  thus  the  Helmholtz  free  energy  is 
strictly  a  function  of  entropy,  allowing  the  nominal  stress  f*  to  be  then  given  by 


r  =  -T 


vkT  rQ 
3 


G'{r0a )  —  a  3^2G'(rQa  1^2) 


(1) 


where  G(r)  =  lnP(r)  and  G'(r)  =  ^ .  The  corresponding  modulus  [/*]  is  computed  from  the 

relation, 


=  r 

a  —  a~ 2 

For  small  strains  (a  — >•  1),  the  modulus  in  (2)  approaches  Young’s  modulus.  Based  on  the  previous 
discussion,  an  expression  for  P(r)  which  is  both  accurate  and  easily  manipulated  mathematically 
is  essential.  The  significance  of  these  implications  as  it  relates  to  macroscopic  stiffness  predictions 
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is  discussed  in  [22], 

In  addition  to  developing  an  appropriate  simulation  model  for  Nafion  polymer  chain  conforma¬ 
tions  as  it  relates  to  rotational  isomeric  state  (RIS)  theory,  another  significant  goal  of  this  paper 
is  to  investigate  and  compare  several  different  methodologies  for  estimating  the  probability  density 
function  P(r).  In  Mark  and  Curro  [12],  a  cubic  spline  approach  is  used  to  estimate  P(r).  Here, 
we  present  two  alternative  methods  for  estimating  P(r)  using  the  Johnson  translation  system  of 
distributions  [10]  and  Bezier  distributions  [21].  While  Bezier  curves  offer  a  very  flexible  approach  to 
density  estimation  with  an  open-ended  parameterization  that  allows  for  the  estimation  of  density 
functions  with  multiple  modes,  our  results  show  that  the  estimates  of  P(r)  obtained  by  fitting  an 
appropriate  Johnson  distribution  to  the  data  are  more  intuitive  than  those  obtained  using  Bezier 
and  cubic  spline  approaches  for  the  following  reasons.  First,  it  is  possible  to  write  down  an  ex¬ 
plicit  functional  form  for  the  Johnson  density  function  that  is  simple  to  differentiate.  This  is  a 
crucial  property  since  the  function  will  be  used  as  an  input  to  the  mathematical  model  (1)  and 
(2)  to  compute  a  Young’s  modulus.  The  functional  form  of  a  Bezier  density  is  considerably  more 
complicated  than  the  Johnson  and  whereas  a  cubic  spline  density  function  has  a  simple  form,  its 
piecewise  definition  results  in  considerably  more  parameters  to  estimate  than  the  globally-defined 
Johnson  density  function.  Second  when  fitting  a  Johnson  density  function  to  the  simulated  data, 
one  continuous  curve  is  fit  to  the  data  so  that  infinite  levels  of  differentiability  are  guaranteed. 
When  fitting  a  cubic  spline  density,  however,  a  piecewise-continuous  curve  is  fit  to  the  data  —  a 
property  that  may  produce  alternating  second  derivative  signs,  which  ultimately  can  lead  to  neg¬ 
ative  stiffness  predictions,  as  shown  in  [22].  Furthermore,  the  piecewise  nature  of  the  cubic  spline 
fit  makes  it  difficult  to  accurately  estimate  the  tails  of  the  distribution. 

Finally,  the  parameters  of  the  Johnson  density  have  a  specific  meaning  (that  is,  they  are  either 
shape,  location,  or  scale  parameters).  As  a  result,  it  may  be  possible  to  find  a  relationship  between 
the  parameters  of  a  Johnson  density  function  and  material  stiffness.  In  [22],  the  results  of  a 
sensitivity  analysis  for  the  parameters  of  the  Johnson  distribution  and  the  corresponding  effect  on 
material  stiffness  is  presented.  Establishing  a  consistent  relationship  between  these  parameters  and 
the  corresponding  stiffness  predictions  would  first  serve  to  extend  RIS  theory  to  stiffness  predictions, 
and  may  ultimately  also  serve  as  a  first  step  toward  the  custom  design  of  material  stiffness  at  the 
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synthesis  stage. 

The  rest  of  this  paper  is  organized  as  follows.  The  development  of  the  simulation  model  for  the 
atomic  structure  of  the  Nation  polymer  backbone  is  described  in  Section  2,  including  details  of  the 
geometry  of  the  Nation  backbone  chain  and  the  rotational  isomeric  scheme  used  to  place  backbone 
atoms  in  three-dimensional  space.  Section  3  includes  descriptions  of  the  cubic  spline,  Johnson,  and 
Bezier  methods  for  density  estimation  as  well  as  the  results  of  applying  these  three  methods  to  the 
simulation-generated  r  values. 

2.  The  Model 

This  paper  focuses  on  the  development  of  a  Monte  Carlo  simulation  model  for  the  atomic  structure 
of  the  Nation  polymer  backbone  chain  with  the  ultimate  goal  of  finding  the  probability  distribution 
of  side  chain  cluster  interaction  distances.  Each  iteration  of  the  simulation  produces  a  single  Nation 
molecule  backbone  within  a  cubic  grid  of  dimension  (5000  A)3.  The  starting  point  of  each  chain  is 
randomly  specified  and  the  backbone  is  dynamically  constructed  according  to  a  commonly  assumed 
geometry  and  under  the  restrictions  that  it  cannot  exceed  the  bounds  of  the  initial  grid  or  collide 
with  the  embedded  clusters. 

2.1.  Geometry  of  the  Nation  Backbone  Chain 

Bulk  Nation  material  is  a  complicated  network  of  molecules.  As  a  polymer,  an  individual  Nation 
molecule  is  composed  of  the  monomer  shown  in  Figure  2,  which  is  repeated  m  times.  The  — (CF2  — 
CF)(CF2  —  CF2)n-  portion  of  the  formula  is  collectively  referred  to  as  the  backbone  of  the  polymer. 
The  -(CF2  —  CF2)n-  portion,  or  Polytetrafluoroethylene  (PTFE),  is  commonly  known  as  Teflon. 
The  portion  of  the  formula  branching  off  at  the  second  carbon  atom  of  the  backbone  monomer  is 
referred  to  as  a  side  chain.  The  SO3  ionic  group  is  hydrophilic  and  generally  interacts  with  the 
nearest  cluster. 
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-  (CF2-CF)-(CF2-CF2)n  - 


o-cf2-cf-o-cf2-cf2-so3- 

cf3 

Figure  2:  Structural  formula  of  Nafion  monomer. 

2.1.1.  Ranges  for  Nafion  Monomer  Variables 

Since  the  structural  formula  for  Nafion  has  not  been  published,  the  values  of  n  and  rn  are  sampled 
from  a  known  range.  The  typical  Nafion  composition  is  87/13,  meaning  that  for  every  87  (CF2CF2) 
groups,  there  are  13  (CF2CF)  groups  comprising  the  polymer  backbone.  The  actual  values  of  n 
used  in  the  simulation  are  drawn  from  the  discrete  probability  distribution  on  the  range  [5,. .  .,11] 
as  summarized  in  Table  1.  The  probabilities  in  Table  1  were  selected  so  that  the  distribution  for  n 
would  be  nearly  symmetric  with  a  mean  of  approximately  87/13  ~  7. 

Preliminary  research  also  suggests  that  there  exists  substantial  variation  in  chain  length  for 
Nafion  polymers.  The  total  molecular  weight  bounding  cases  of  150,000  and  250,000  g/rnol  were  em¬ 
ployed  to  arrive  at  a  reasonable  range  for  rn  to  account  for  this  variation  of  molecular  chain  length. 
Using  the  standard  atomic  masses  for  carbon  (12.0107),  fluorine  (18.9984),  oxygen  (15.9994),  and 
sulfur  (32.0650),  the  mass  of  each  (CF2CF)  portion  of  the  polymer  backbone  with  the  attached 
side  chain  is  443.1161  g/rnol.  When  taken  in  conjunction  with  the  87/13  composition  of  Nafion, 
approximately  40%  of  the  total  molecular  weight  is  produced  by  the  (CF2CF)  portions  of  the  poly¬ 
mer  backbone  along  with  the  attached  side  chains.  Upon  examination  of  the  bounding  cases,  we 


Table  1:  Discrete  probability  distribution  for  n  values. 


n 

Probability 

5 

0.1 

6 

0.2 

7 

0.3 

8 

0.2 

9 

0.1 

10 

0.05 

11 

0.05 
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see  that  0.40x150,000/443.1161  «  135  and  0.40x250,000/443.1161  «  225.  Hence,  the  m  values 
used  in  the  simulation  were  uniformly  drawn  from  integers  in  the  range  [135,  . .  .,225]. 

While  this  is  the  target  range  for  the  repeated  monomer  units,  the  simulated  chains  are  termi¬ 
nated  at  much  shorter  lengths  due  to  chain  collision  with  a  cluster.  Actual  simulated  chain  lengths 
m  ranged  between  1  and  45  with  a  typical  value  of  4.  Whereas  early  termination  clearly  leads  to 
non-physical  predictions  of  total  chain  lengths,  it  is  not  believed  that  early  termination  presents  a 
significant  problem  in  the  estimation  of  P{r)  because  it  is  the  distance  between  cluster  interactions 
which  dictate  the  elastic  properties  of  the  material. 

2.1.2.  Bond  Length  and  Angles  of  Bond  Orientation 

The  initial  model  focuses  on  characterizing  the  backbone  of  the  polymer.  From  standard  chemistry 
principles,  the  length  of  a  carbon-carbon  bond  from  the  centers  of  both  atoms  is  known  to  be 
l  =  1.53  A  [8].  This  fixed  bond  length  was  used  throughout  the  simulation  as  the  distance  between 
backbone  chain  atoms. 

Since  the  simulation  is  performed  in  three-dimensional  space,  there  are  two  angles  required  for 
the  placement  of  any  backbone  atom.  These  are  the  valence  angle  9  between  successive  bonds  and 
the  angle  of  rotation  (f>  for  each  bond.  For  our  purposes,  the  valence  angle  (like  the  bond  length)  is 
considered  to  be  a  fixed  quantity  for  all  pairs  of  successive  bonds.  Due  to  the  large  proportion  of 
PTFE  in  the  backbone,  well-established  PTFE  values  are  applied  throughout.  In  particular,  the 
valence  angle  9  is  fixed  at  116°.  The  parameter  of  experimental  interest  is  the  angle  of  rotation  <\>. 

As  depicted  in  Figure  3,  the  angle  of  rotation  cj)j  for  bond  j  is  measured  in  relationship  to  the 
plane  formed  by  the  three  preceding  carbon  atoms  of  the  backbone  chain.  Given  this  plane  and 
the  fixed  valence  angle  9,  an  atom  i  may  be  placed  anywhere  along  the  base  of  the  cone  formed  by 
the  rotation  of  the  jth  bond.  To  be  described  in  further  detail  in  Section  2.2,  Teflon  is  assumed  to 
have  four  possible  values  of  ipj'.  ±15°  and  ±120°.  Hence,  an  atom  i  may  be  placed  in  any  of  four 
possible  positions  along  the  base  of  the  cone  shown  in  Figure  3.  Specifically,  <pj  is  measured  in  a 
clockwise  manner  from  the  point  on  the  base  of  the  cone  lying  in  the  plane  formed  by  the  three 
preceding  atoms  when  in  the  fully  extended  planar  zigzag  configuration,  which  would  have  a  value 
of  fa  =  0°  [7]. 
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wi-3 

Figure  3:  Spatial  geometry  of  a  Nation  molecule. 

For  simplicity,  when  placing  an  atom  i,  the  simulation  temporarily  performs  the  necessary 
rotations  and  translations  to  place  the  backbone  atom  i  —  1  at  the  origin,  atom  i  —  2  along  the 
negative  x-axis  of  the  Cartesian  coordinate  grid,  and  atom  i  —  3  in  the  negative  x-negative  y 
quadrant  of  the  x-y  plane.  Once  in  this  configuration,  atom  i  may  be  placed  locally  by  randomly 
selecting  an  angle  cf>j  and  measuring  clockwise  from  the  point  on  the  base  of  the  cone  in  the  positive 
x-positive  y  plane,  as  shown  in  Figure  4.  Specifically,  the  local  Cartesian  coordinates  of  atom  i  are 
defined  by 

Xj  £  ■  cos  (O') 

yi  =  l  ■  cos (—<t>j)  ■  sin(0/) 

Zi  £  ■  sin (—<j>j)  •  sin(0') 

where  £  is  the  carbon-carbon  bond  length  of  1.53  A  and  0'  is  the  supplementary  angle  to  9  having 
the  value  of  180  —  0  =  64°.  After  placing  atom  i  locally,  the  inverse  of  the  rotations  and  translations 
applied  to  the  positions  of  atoms  i  —  1,  i  —  2,  and  i  —  3  to  achieve  the  local  configuration  in  Figure 
4  are  performed  on  atom  i  so  that  it  is  in  the  correct  position  with  respect  to  the  remainder  of  the 
backbone  chain  before  temporarily  rotating  atoms  i  —  1,  i  —  2,  and  i  —  3. 

2.2.  Rotational  Isomeric  Scheme 

Although  the  conformational  characteristics  of  PTFE  or  Teflon  (the  predominant  structure  in  the 
Nation  backbone)  are  well  represented  by  a  three-state  model  with  rotation  angles  of  one  trans 
(4>  =  0°)  and  two  gauche  ((f)  =  ±115°)  conformations,  it  has  been  shown  that  a  four-state  model  is 
more  physically  accurate  [3].  Therefore,  this  simulation  adopts  the  four-state  model  with  rotation 
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Figure  4:  Local  geometry  of  a  Nation  molecule  after  temporarily  rotating  atoms  i  —  1,  i  —  2,  and 
i  —  3  for  the  purpose  of  placing  atom  i. 

angles  of  two  trans  (<f>  =  ±15°)  and  two  gauche  (cf)  =  ±120°)  conformations. 

2.2.1.  Statistical  Weight  Matrices 

From  [7],  it  is  known  that  the  given  rotational  state  of  any  bond  in  a  polymer  depends  on  the 
rotational  state  of  its  neighboring  bonds  .  The  first-neighbor  dependence  was  incorporated  into 
this  simulation  by  invoking  a  statistical  weight  parameter  u(a,rj,j )  =  where  a  is  the  state  of 

the  preceding  bond,  g  is  the  state  of  the  current  bond,  and  j  is  the  location  of  the  bond  within  the 
molecular  chain.  Therefore,  adapting  from  the  values  for  the  statistical  weight  parameters  arrived 
at  by  Bates  and  Stockmayer  in  [3],  we  employ  the  matrix  Uj  =  [ua,ri,j\  for  bond  j.  where 
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The  values  of  a  =  0.2,  a'  =  2.0,  uj  =  0.2,  and  (3  =  0.5  were  identified  by  comparing  the  mean  of 
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the  range  for  each  of  the  values  provided  in  [3]  and  the  associated  energy  values.  Note  that  all  of 
the  above  matrices  are  indexed  in  the  order  t+  (15°),  g+  (120°),  g~  (-120°),  and  t~  (-15°)  on  both 
the  rows  and  columns.  Furthermore,  since  the  first  bond  is  placed  randomly,  we  have  included  U2 
to  place  the  second  bond  of  the  backbone,  U3  to  place  the  third  bond  of  the  backbone,  and  Un  to 
place  the  last  bond  of  the  entire  backbone  chain.  Hence  Uk  is  used  to  place  bonds  4  <  k  <  N  —  1, 
where  N  =  2 (m  +  np)  —  1  is  the  total  number  of  bonds  in  the  backbone  chain  and  np  is  the 
pth  randomly  generated  value  of  n  from  the  range  [5, . . . ,  11],  as  described  in  Section  2.1.1. 

2.2.2.  Conditional  Probability  Matrices 

The  stochastic  element  of  the  Monte  Carlo  simulation  of  the  atomic  structure  of  Nation  focuses 
on  determining  the  angle  of  rotation  (j)  for  each  bond.  However,  the  statistical  weight  matrices 
Uj,  2  <  j  <  N  derived  in  the  previous  section  are  non- normalized.  Therefore,  we  must  use  the 
similarity  transformation 

Qj  =  £j.\  UjDj_\, 

where  £j  1  is  the  largest  eigenvalue  of  Uj  and  Dj\  is  the  diagonal  matrix  of  the  elements  of  the 
eigenvector  associated  with  £j  1,  to  normalize  each  of  the  Uj  to  obtain  the  conditional  probability 
matrices  [8], 
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As  in  the  statistical  weight  matrices  of  the  previous  section,  j  represents  the  location  of  the 
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bond  in  the  entire  molecular  chain  and  each  Qj  is  indexed  in  the  order  t+  (15°),  g+  (120°),  g~ 
(-120°),  and  t~  (-15°)  on  both  the  rows  and  columns.  Note  that  the  rows  of  the  above  matrices  do 
not  sum  exactly  to  unity  since  each  value  has  been  truncated  in  the  interest  of  space.  It  is  these 
conditional  probability  matrices  which  are  employed  to  choose  the  angle  of  rotation  <f>  for  all  bonds 
in  the  simulation  based  on  the  rotational  state  of  the  preceding  bond,  where  in  the  actual  code 
the  untruncated  values  are  used  so  that  each  row  sums  exactly  to  unity.  For  example,  suppose  we 
are  attempting  to  place  the  last  bond  of  a  backbone  chain  where  the  preceding  bond  is  in  the  t+ 
state.  This  indicates  that  the  first  row  of  Qn  contains  the  appropriate  conditional  probabilities  to 
select  the  particular  cj)  value  for  the  current  bond.  A  random  number  generator  is  invoked  to  select 
a  random  number  between  0  and  1.  Assume  a  value  of  0.4  is  returned.  Since  0.278  <  0.4  <  (0.278 
+  0.482),  we  have  randomly  selected  the  value  =  120°. 

2.3.  Initial  Grid  Orientation 

As  noted  previously,  the  space  within  which  the  molecular  chains  are  constructed  is  of  dimension 
(5000  A)3.  The  standard  Cartesian  coordinate  system  is  oriented  so  that  the  range  of  values  along 
each  of  the  x,  y,  and  z  axes  is  [-2500,  2500].  These  dimensions  were  chosen  so  that  the  longest 
possible  molecule,  according  to  the  maximum  m  and  mean  n  values  of  the  Nafion  formula,  when 
fully  extended,  could  be  contained  within  the  grid. 

It  is  of  particular  interest  in  this  investigation  to  examine  the  effects  of  embedded  inclusions  upon 
the  cluster  communication  distances.  In  this  simulation,  the  embedded  inclusions  are  hydrophilic 
clusters  in  which  the  side  chains  tend  to  terminate  and  which  are  assumed  to  already  exist  as 
though  the  grid  was  populated  with  other  Nafion  molecules.  This  simplified  model  considers  the 
clusters  to  be  spherical  in  shape  and  uniformly  distributed  throughout  the  grid  [9].  Since  it  is 
difficult  to  arrive  at  a  truly  uniform  distribution  of  spheres  on  a  rectangular  grid  where  the  centers 
of  the  spheres  are  equidistant  from  each  other,  two  possible  grid  orientations  were  considered: 
rectangular  (analogous  to  cubic  crystalline  structure)  and  staggered  (analogous  to  hexagonal  close 
packed  (HCP)  crystalline  structure),  as  shown  in  Figure  5. 

The  center-to-center  distance  between  clusters  in  either  orientation  is  determined  by  the  choice 
of  cluster  radius,  which  was  considered  to  be  identical  for  all  clusters  populating  the  grid,  as  well  as 
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(a) 


(b) 


Figure  5:  (a)  Rectangular  grid  orientation,  and  (b)  staggered  grid  orientation. 


the  choice  of  volume  fraction  of  the  entire  grid  occupied  by  the  clusters.  The  radius  of  the  spherical 
clusters  and  the  volume  fraction  of  the  total  grid  they  occupy  were  varied  to  examine  several 
specific  cases.  The  first  of  these  is  the  case  where  the  1100  Equivalent  Weight  Nafion  material  is 
fully  hydrated  with  water  and  sodium  is  used  as  the  cation.  Under  these  conditions,  the  clusters 
each  have  a  radius  of  21  A  and  the  total  volume  fraction  of  the  clusters  is  30%,  based  on  results 
from  Bar-Cohen  [2].  From  these  characteristics,  the  center-to-center  distance  between  clusters  for 
both  the  rectangular  and  staggered  orientations  is  2.4xradius  =  50.4  A.  In  the  second  case,  we 
assumed  that  the  1100  Equivalent  Weight  Nafion  material  is  fully  hydrated  with  water  and  lithium 
is  used  as  the  cation.  These  conditions  result  in  a  radius  of  23  A  and  the  total  volume  fraction  of 
the  clusters  to  be  38%,  based  on  results  from  Bar-Cohen  [2].  These  values  imply  a  center-to-center 
distance  between  clusters  for  both  the  orientations  to  be  2.2xradius  =  50.6  A. 

2.4.  Variations  for  Backbone  Atoms  with  Side  Chains 

2.4.1.  Atom  Placement 

In  Section  2.2.2,  a  method  for  selecting  the  angle  of  rotation  <f>  for  each  bond  was  discussed.  For 
those  atoms  with  attached  side  chains,  however,  we  propose  the  following  alternative  method  for 
selecting  cf).  Due  to  the  hydrophilic  nature  of  the  side  chain  terminal  groups  in  an  otherwise 
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hydrophobic  molecule,  there  exists  a  tendency  for  the  terminal  groups  and  the  water  absorbed  by 
the  material  to  cluster  together.  We  choose  to  incorporate  this  phenomenon  into  the  simulation  by 
overriding  the  conditional  probability  matrices  used  to  choose  0  (as  described  in  Section  2.2.2)  and 
instead  use  the  0-value  that  would  place  the  backbone  atom  with  the  attached  side  chain  closest  to 
an  embedded  cluster.  This  is  accomplished  by  forming  all  possible  orientations  produced  by  each 
of  the  ±15°  and  ±120°  rotation  angles  and  computing  the  distances  to  the  nearest  clusters.  Once 
the  minimum  distance  is  found,  the  associated  angle  of  rotation  0  is  used  to  actually  place  the 
backbone  atom  with  the  attached  side  chain  using  the  same  method  described  in  Section  2.1.2. 

2.4.2.  Achieving  Side  Chain  Cluster  Communication 

It  was  noted  previously  that  the  goal  in  this  characterization  technique  is  to  quantify  the  distance 
between  cluster  interaction  points  along  the  Nafion  molecular  backbone.  At  every  backbone  atom 
with  an  attached  side  chain,  there  is  an  opportunity  for  a  cluster  interaction  point,  meaning  that 
the  attached  side  chain  terminal  group  comes  into  communication  with  an  embedded  cluster.  This 
cluster  communication  is  assumed  to  occur  when  the  backbone  atom  with  an  attached  side  chain 
is  within  a  distance  of  8  A  of  an  embedded  cluster.  The  cluster  communication  distance  of  8  A  was 
determined  by  examining  the  composition  of  the  attached  side  chains.  Although  they  may  take 
on  a  number  of  conformations  and  hence  lengths,  the  maximum  length  is  approximated  as  8  A. 
Therefore,  the  backbone  atoms  which  achieve  side  chain-cluster  communication  are  tracked  such 
that  successful  cluster  communication  is  assumed  to  act  as  a  crosslink.  Our  end  result  is  comprised 
of  computing  the  distances  between  these  subsequent  backbone  atoms,  as  well  as  the  endpoints  of 
the  molecular  chain. 

Quantifying  the  rate  of  successful  cluster  interactions  —  that  is,  the  ratio  of  successful  cluster 
interactions  to  the  total  number  of  possible  interactions  —  is  of  scientific  interest.  Table  2  summa¬ 
rizes  the  successful  cluster  interaction  rates  when  a  cluster  communication  distance  of  8  A  is  used. 
From  this  table  we  see  that  the  rates  produced  by  the  four  different  cases  agree  with  simple  intu¬ 
ition.  Since  the  embedded  clusters  of  the  staggered  orientation  tend  to  be  slightly  closer  together 
than  the  rectangular  orientation,  the  result  is  a  higher  interaction  success  rate.  Similarly,  since 
the  case  where  Nafion  is  fully  hydrated  with  lithium  as  the  cation  has  a  larger  volume  fraction  of 
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Table  2:  Successful  cluster  communication  rates. 


Grid  Orientation 

Cation 

Rate 

Rectangular 

Sodium 

75% 

Rectangular 

Lithium 

83% 

Staggered 

Sodium 

86% 

Staggered 

Lithium 

92% 

clusters,  this  would  also  suggest  higher  success  rates  due  to  tighter  cluster  packing.  Furthermore, 
although  we  have  employed  the  cluster  communication  distance  of  8  A,  only  approximately  10%  of 
the  successful  cluster  interactions  occur  when  the  backbone  atom  with  the  attached  side  chain  is 
between  6.5  A  and  8  A  from  the  nearest  cluster.  Hence  the  bulk  of  the  simulated  cluster  interac¬ 
tions  occur  at  minimal  distances,  which  is  critical  due  to  the  fact  that  8  A  is  an  assumed  maximum 
value. 


3.  Simulation  Results 

The  simulation  model  was  run  for  each  of  the  following  four  cases:  sodium  as  the  cation  with 
a  rectangular  grid  orientation,  sodium  as  the  cation  with  a  staggered  grid  orientation,  lithium 
as  the  cation  with  a  rectangular  grid  orientation,  and  lithium  as  the  cation  with  a  staggered  grid 
orientation.  For  each  of  the  four  cases,  a  data  set  of  approximately  10,000  r  values  was  generated  by 
simulating  possible  conformations  of  the  Nafion  backbone,  using  the  method  described  in  Section  2. 
We  then  estimated  the  probability  density  function  P(r )  using  two  methods:  Johnson  distributions 
and  Bezier  functions.  In  the  Sections  3.2  and  3.3,  we  will  detail  the  steps  followed  to  estimate  P{r) 
using  each  of  these  two  methods.  First,  however,  we  will  describe  the  traditionally  applied  cubic 
spline  method  for  density  estimation. 

3.1.  Fitting  Data  with  Cubic  Splines 

In  their  work  on  short  chain  polymers  displaying  rubber  elasticity,  Mark  and  Curro  [12]  use  cubic 
splines  to  estimate  the  probability  density  function  P(r)  for  chain  lengths  of  polyethylene  and 
polydimethylsiloxane.  The  simulation-generated  r  values  are  divided  into  groups,  or  bins,  and 
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a  knot  value  is  specified  for  each  bin  so  that  on  each  interval  (delimited  by  the  knots),  a  cubic 
polynomial  is  fit  to  the  data.  If  we  define  the  indicator  function  Lj(r)  for  the  subinterval  (ki,  fcj+i], 
where  ki  is  the  ith  knot,  as 


U(r) 


1,  for  ki  <r  <  fcj+i, 
0,  otherwise, 


(3) 


then  the  result  is  a  piecewise-polynomial  estimate  of  P(r)  of  the  form 

p(r)  = 

i=  1 

where  n*  is  the  number  of  knots,  A,  B,  C,  and  D  are  the  estimated  coefficients,  and  K  is  a  normal¬ 
izing  constant  so  that  the  function  P(r)  given  in  (4)  integrates  to  unity. 

The  advantage  to  using  cubic  splines  for  density  estimation  is  that  the  resulting  function  has 
a  simple  form  and  can  easily  be  applied  in  mathematical  models  to  describe  various  material 
properties.  One  challenge  with  this  approach  is  that  cubic  polynomials  are  often  not  sufficiently 
flexible  to  accurately  model  the  tail  behavior  of  the  density.  Additionally,  the  shape  of  the  estimated 
density  is  highly  dependent  on  the  number  and  location  of  the  knots.  For  instance,  if  the  number 
of  knots  is  large,  then  the  corresponding  estimate  of  P(r)  may  be  excessively  noisy.  Finally,  the 
details  of  the  fitting  methodology  are  typically  not  well-defined  and  often  a  trial-and-error  approach 
is  necessary  to  determine  the  number  and  location  of  the  knots. 

In  Figure  6,  a  plot  of  the  estimated  cubic  spline  density  function  for  the  sodium/rectangular 
case  is  shown  with  a  histogram  of  the  10,000  simulated  r  values  superimposed.  From  this  plot 
it  is  clear  that  it  is  possible  to  achieve  an  excellent  estimate  of  the  distribution  of  chain  length 
using  cubic  splines.  However,  as  indicated  in  [22],  minor  variation  in  the  number  and  placement 
of  knots  can  have  a  non-negligible  impact  on  the  predicted  stiffness  and,  furthermore,  the  cubic 
spline  estimates  of  P(r)  often  lead  to  non-physical  (negative)  stiffness  predictions.  We  believe  the 
negative  stiffness  predictions  are  partially  the  result  of  the  piecewise  nature  of  the  estimated  cubic 
spline  density  function,  which  can  lead  to  unpredictable  behavior  in  the  second  derivative  of  P(r). 
In  the  following  subsections,  we  propose  two  alternative  methods  for  estimating  the  distribution  of 


Ai(r  —  A:*)3  +  Bi(r  —  ki )2  +  Ci(r  —  ki)  +  D{ 


(4) 
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Figure  6:  Estimated  cubic  spline  PDF  for  the  sodium/rectangular  case. 


chain  lengths. 

3.2.  Fitting  Data  with  Johnson  Distributions 

The  Johnson  translation  system  of  distributions  [10]  is  comprised  of  four  flexible  families  of  dis¬ 
tributions  including  lognormal,  unbounded,  bounded,  and  normal  which  cover  a  wide  variety  of 
distribution  shapes.  It  is  this  flexibility  which  makes  the  Johnson  system  an  attractive  option  to 
estimate  P(r)  from  the  simulated  data  for  polymer  chain  length. 

If  X  is  the  continuous  random  variable  whose  distribution  is  to  be  estimated,  the  cumulative 
distribution  function  (CDF)  of  X  is  F(x)  =  Pr{X  <  x}  and  the  probability  density  function  (PDF) 
of  X  is  P{x)  =  F'{x).  Johnson  developed  a  method  for  estimating  the  unknown  density  P{x )  by 
transforming  the  random  variable  X  into  a  standard  normal  random  variable, 


Z  =  7  +  Sf 


X-tj 

A 


where  Z  is  a  standard  normal  random  variable,  7  and  5  are  shape  parameters,  A  is  a  scale  parameter, 
£  is  a  location  parameter,  and  /  is  one  of  the  following: 
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f(v)  =  { 


(5) 


My), 

ln(y  +  \Jy 2  +  1), 

Mi**), 

,  y> 


for  the  Sl  (lognormal)  family, 
for  the  Su  (unbounded)  family, 
for  the  Sb  (bounded)  family, 
for  the  S]\r  (normal)  family. 


The  corresponding  PDF  for  X  is 


p(x )  = 


(6) 


for  all  x  €  H  where  /'(•)  is  the  first  derivative  of  the  function  /(•)  in  (5)  and  the  (closed)  support 
of  the  distribution  is 


[£,+oo) 

(— oo,+oo) 

[£,  £  +  M 

(-oo,+oo) 


for  the  Sl  (lognormal)  family, 
for  the  Sjj  (unbounded)  family, 
for  the  Sb  (bounded)  family, 
for  the  S]\r  (normal)  family. 


The  software  package  FITTR1  [24]  incorporates  several  methods  for  fitting  Johnson  distribu¬ 
tions,  including  moment  matching,  percentile  matching,  ordinary  least  squares  (OLS),  diagonally 
weighted  least  squares  (DWLS),  Li-norm  and  Loo-norm  estimation  for  each  of  the  four  families. 
Descriptive  statistics  (including  mean,  standard  deviation,  skewness,  and  kurtosis)  of  the  sample 
data,  goodness-of-fit  information,  as  well  as  tables  of  the  empirical  and  fitted  CDF  and  PDF  for 
plotting  are  available  outputs  of  the  FITTR1  package. 

Among  the  goodness-of-fit  information  provided  as  output  by  FITTR1  is  the  Kolmogorov- 
Smirnov  statistic,  which  indicates  the  maximum  discrepancy  between  the  empirical  and  fitted 
distributions  as  a  percentage  value.  Note  that  in  this  context,  empirical  distribution  refers  to 
the  simulated  distribution,  as  true  experimental  data  is  presently  unattainable.  The  Kolmogorov- 
Smirnov  statistic  was  used  as  a  metric  to  determine  the  accuracy  of  the  fitted  distributions.  For 
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each  of  the  four  simulated  cases,  estimates  of  P(r)  were  generated  using  various  combinations 
of  Johnson  families  (primarily  the  bounded  and  unbounded)  and  fitting  methods.  Based  on  the 
Kolmogorov-Smirnov  statistic  and  plots  of  the  estimated  Johnson  PDF  and  CDF  versus  the  data, 
one  fit  was  selected  as  the  “best”  fit  for  each  of  the  four  cases  (see  Figures  7-10).  Table  3  provides 
the  corresponding  parameter  estimates  and  the  Kolmogorov-Smirnov  goodness-of-fit  statistic  (K-S 
g.o.f.)  for  the  Johnson  distributions  shown  in  Figures  7-10. 

Although  our  goal  is  to  determine  an  estimate  of  P(r),  the  PDF  for  the  cluster  interaction 
distances  r,  we  select  the  best  fit  primarily  by  comparing  plots  of  the  empirical  CDF  versus  the 
fitted  Johnson  CDF  since  the  empirical  PDF,  which  is  often  represented  as  a  histogram,  is  highly 
subjective  and  the  shape  can  vary  greatly  depending  on  the  choice  of  histogram  cell  width.  For 
example,  if  the  width  of  the  histogram  bars  are  taken  to  be  very  small,  the  resulting  estimate  of 
the  PDF  may  be  excessively  noisy  or  bumpy  (that  is,  it  may  contain  extraneous  bumps  or  modes 
that  are  a  function  of  the  randomness  of  the  simulation  and  may  not  actually  be  features  of  the 
true  underlying  PDF).  Conversely,  if  the  width  of  the  histogram  bars  are  taken  to  be  very  large, 
the  resulting  estimate  of  the  PDF  may  be  over  smoothed  and  important  features  of  the  true  PDF 
may  be  lost. 

When  studying  the  plots  of  the  empirical  versus  the  fitted  CDF,  we  consider  two  main  regions: 
the  middle  section  (10  <  r  <  30)  and  the  tail  sections  (r  <  10  and  r  >  30).  Of  greater  importance 
is  the  closeness  of  fit  in  the  middle  section  since  this  region  characterizes  the  bulk  of  the  cluster 
communication  distances.  As  seen  in  Figures  7-10,  we  were  able  to  successfully  fit  the  middle 
portion  of  the  data  using  a  Johnson  unbounded  distribution  for  each  of  the  four  cases.  Accuracy 
while  fitting  the  tail  sections  is  also  desired,  however  it  is  difficult  in  this  case  to  get  a  close  fit 
in  the  lower  tail  since  the  fitted  Johnson  density  is  continuous  and  the  data  in  the  lower  tail 


Table  3:  Parameter  estimates  for  the  fitted  Johnson  unbounded  distribution  for  each  of  the  four 
simulated  cases. 


Cation 

Grid  Orientation 

Family 

Fitting  Method 

Pr(r  <  0) 

K-S  g.o.f. 

7 

8 

A 

£ 

Sodium 

Rectangular 

Su 

DWLS 

0.0126 

0.0471 

-0.4963 

0.9736 

5.301 

15.44 

Sodium 

Staggered 

Su 

OLS 

0.0356 

0.0415 

-0.0247 

0.7758 

3.562 

17.46 

Lithium 

Rectangular 

Su 

Li-norm 

0.0363 

0.0420 

-0.0986 

0.7463 

3.596 

17.28 

Lithium 

Staggered 

Su 

OLS 

0.0463 

0.0543 

0.3108 

0.9147 

4.229 

18.44 
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is  discrete.  For  those  r  values  computed  for  backbone  subchains  with  a  very  small  number  of 
bonds  (less  than  5),  there  is  only  a  small  number  of  combinations  of  atom  placement  since  there 
are  only  four  angles  of  rotation  to  choose  from.  For  example,  if  the  formation  of  the  backbone 
chain  terminates  after  the  placement  of  the  second  carbon  atom,  then  a  value  of  r  =  1.53  will  be 
recorded,  corresponding  to  the  length  of  one  carbon-carbon  bond.  Similarly,  if  the  backbone  chain 
terminates  after  the  placement  of  the  third  carbon  atom,  then  a  value  of  r  =  2.595,  corresponding 
to  the  distance  between  the  ends  of  2  carbon-carbon  bonds,  will  be  recorded.  If  the  backbone  chain 
terminates  after  the  placement  of  the  fourth  atom,  then  there  are  only  four  possible  resulting  r 
values  corresponding  to  each  possible  value  of  the  rotation  angle  cj).  As  the  number  of  bonds  in 
a  backbone  subchain  increases,  the  distribution  of  r  values  begins  to  look  more  continuous  since 
there  are  increasingly  more  combinations  of  atom  placement  and  hence  chain  lengths. 


(a)  (b) 

Figure  7:  Estimated  PDF  (a),  and  CDF  (b)  for  the  so dium/rect angular  case  using  a  Johnson 
unbounded  density  function. 
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(a)  (b) 

Figure  8:  Estimated  PDF  (a),  and  CDF  (b)  for  the  sodium/staggered  case  using  a  Johnson  un¬ 
bounded  density  function. 


(a)  (b) 


Figure  9:  Estimated  PDF  (a),  and  CDF  (b)  for  the  lithium/rectangular  case 
unbounded  density  function. 
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(a)  (b) 

Figure  10:  Estimated  PDF  (a),  and  CDF  (b)  for  the  lithium/staggered  case  using  a  Johnson 
unbounded  density  function. 


21 


Another  critical  drawback  to  the  Johnson  unbounded  fits  in  Figures  7-10  is  that  there  is  a 
significant  nonzero  probability  that  the  fitted  Johnson  density  could  generate  an  r  value  less  than 
zero,  as  seen  in  Table  3.  To  correct  this,  we  attempted  to  fit  a  Johnson  bounded  distribution  to  the 
data  for  each  of  the  four  cases.  The  fitted  CDF  and  PDF  for  the  Johnson  bounded  fits  are  shown 
in  Figures  11-14  and  Table  4  gives  the  K-S  g.o.f.  statistic  and  the  parameter  values  for  each  of  the 
four  cases.  The  lower  bound  for  each  case  was  set  to  zero.  Clearly,  while  the  Johnson  bounded  fits 
have  a  zero  probability  of  generating  a  negative  r  value,  the  fits  are  not  as  accurate  as  the  Johnson 
unbounded  fits.  Since  the  lower  bound  is  fixed,  the  Johnson  bounded  distribution  only  has  three 
parameters  to  estimate,  rendering  it  even  less  flexible  than  the  Johnson  unbounded  distribution. 

Despite  these  drawbacks,  however,  the  Johnson  distribution  still  provides  an  accurate  estimate  of 
the  density  of  polymer  chain  length,  requiring  the  estimation  of  only  four  parameters.  Furthermore, 
the  functional  form  of  the  Johnson  PDF  in  (6)  can  easily  be  used  as  an  input  to  the  macroscopic- 
level  model  given  by  (2)  to  predict  material  stiffness.  As  indicated  in  [22],  the  estimated  Johnson 
density  functions  for  chain  length  result  in  stable  and  repeatable  stiffness  predictions. 


(a) 


(b) 


Figure  11:  Estimated  PDF  (a),  and  CDF  (b)  for  the  sodium/rectangular  case  using  a  Johnson 
bounded  density  function. 
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(a)  (b) 


Figure  12:  Estimated  PDF  (a),  and  CDF  (b)  for  the  sodium/staggered  case  using  a  Johnson 
bounded  density  function. 


(a)  (b) 


Figure  13:  Estimated  PDF  (a),  and  CDF  (b)  for  the  lithium/rectangular  case  using  a  Johnson 
bounded  density  function. 


(a)  (b) 


Figure  14:  Estimated  PDF  (a),  and  CDF  (b)  for  the  lithium/staggered  case  using  a  Johnson 
bounded  density  function. 
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Table  4:  Parameter  estimates  for  the  fitted  Johnson  bounded  (Sb)  distribution  for  each  of  the  four 
simulated  cases. 


Cation 

Grid  Orientation 

Family 

Fitting  Method 

K-S  g.o.f. 

7 

8 

A 

£ 

Sodium 

Rectangular 

SB 

OLS 

0.0797 

2.760 

2.076 

87.96 

0 

Sodium 

Staggered 

Sb 

OLS 

0.1041 

2.721 

2.385 

72.19 

0 

Lithium 

Rectangular 

Sb 

Lj-norm 

0.1050 

3.059 

2.495 

78.23 

0 

Lithium 

Staggered 

Sb 

OLS 

0.1205 

1.910 

2.138 

57.29 

0 

3.3.  Fitting  Data  with  Bezier  Distributions 

Although  the  Johnson  translation  system  of  distributions  is  remarkably  flexible,  it  is  limited  by 
the  fact  that  only  four  parameters  are  available  to  control  the  shape  of  the  fitted  distribution. 
Alternatively,  the  Bezier  distribution  family  offers  an  open-ended  parameterization  technique  which 
enables  it  to  accurately  represent  nearly  any  possible  distributional  shape  [21].  This  characteristic 
enables  the  fitted  Bezier  curve  to  characterize  the  behavior  in  the  tails  of  the  empirical  distribution 
more  precisely. 

Bezier  curves  are  a  special  class  of  spline  curves  used  to  approximate  a  smooth  univariate  func¬ 
tion  on  a  bounded  interval  by  forcing  the  curves  to  pass  through  the  vicinity  of  selected  control 
points  {pfc  =  (xk,  Zk)T  :  k  =  0, 1, . . . ,  n}.  Visually,  the  control  points  can  be  interpreted  as  ‘mag¬ 
nets’,  exerting  a  ‘magnetic  attraction’  which  forces  the  Bezier  curve  towards  them.  In  particular, 
the  first  and  last  control  points  are  interpolated  exactly,  meaning  that  the  Bezier  curve  will  pass 
through  these  points.  It  should  be  noted  that  changing  the  location  of  any  control  point  affects  the 
shape  of  the  entire  curve.  Every  control  point  provides  two  parameters,  specifically  the  Xk  and 
Zk  coordinates.  Hence,  if  a  Bezier  curve  of  degree  n  is  used,  we  have  n  +  1  control  points  and  2 n 
parameters. 

A  Bezier  curve  of  degree  n  with  control  points  {po,  pi, . . . ,  p„}  is  parametrically  written  as 
P  (t)  =  [Px{t]n,x),Pz(t]n,z)}T 

n 

=  ^Bn,k{t) pk  fortG  [0,1], 

k=0 

where  x  =  (xq,  xi, . . . ,  xn)T  and  z  =  (zq,  z±,  . . . ,  zn)T  respectively  denote  the  vectors  of  x-  and 
z-coordinates  of  the  given  control  points,  and  BU)k(t)  are  the  Bernstein  polynomials  defined  as 
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Bn,k(t)  = 


kl(n-k)ltkO-  -  t)n  k>  for  t  €  [0,1]  and  k  =  0, 1, . . . ,  n, 

0,  otherwise. 

If  X  is  a  (continuous)  Bezier  random  variable  with  bounded  support  [x*,x*],  the  cumulative 
distribution  function  (CDF)  Fx[x(t )]  can  be  approximated  to  arbitrary  accuracy  by  a  Bezier  curve 
with  sufficiently  high  degree  n.  Hence,  the  CDF  is  given  parametrically  by 


P(t)  =  [ Px(t-,n,x),Pz(t-n,z)]T 

=  {x(t),  Fx[x(t)]}  for  all  t  G  [0,1] 

where  the  control  points  {po,  pi, . . . ,  pn}  are  placed  so  as  to  adhere  to  the  basic  requirements 
of  a  CDF — that  is,  Fx(- )  is  monotonically  nondecreasing,  Fx{x*)  =  0  and  Fx(x*)  =  1.  Given 
this  representation  of  the  CDF,  the  corresponding  probability  density  function  (PDF)  fx[x(t)\  = 
F'x[x(t)]  is  given  parametrically  by 

P*W  =  [■ P*{t-,n,x),P*(t-,n,x,z)]T  ^ 

=  {x(t),fx[x(t)]}  for  all  t  G  [0,1]. 


Here 


P*(t;n,x,z) 


Pz(t;n  -  1,  A z) 
Px(t;n-  1,  Ax)’ 


Ax  =  (Ax0,  •  •  • ,  Axn  -  1)T, 


Az  =  (Az0, . . . ,  Azn  -  1)T, 


and  Axfc  =  x^+i  —  x^,  and  A z^  =  zj~+ 1  —  Zk{k  =  0,1,..., n  —  1)  are  taken  to  represent  the 
corresponding  differences  of  the  x-  and  z-  coordinates  of  the  control  points  in  the  parametric 
representation  of  the  CDF. 

To  fit  a  Bezier  distribution  to  the  simulated  data  for  each  of  the  four  cases,  we  use  the  graphical 
software  package  PRIME.  Among  the  fitting  methods  available  in  PRIME  are  moment  match¬ 
ing,  percentile  matching,  maximum  likelihood,  least  squares,  Li-  and  Loo-nornr  estimation.  After 
loading  the  data,  the  software  can  estimate  the  number  of  control  points  necessary  to  attain  an  ad- 
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equate  fit,  or  the  user  can  interactively  choose  both  the  number  and  location  of  the  control  points. 
Descriptive  statistics,  goodness-of-fit  information  (including  the  Kolmogorov-Smirnov  statistic),  as 
well  as  plots  of  the  empirical  and  fitted  CDF  and  PDF  are  available  outputs  of  the  PRIME  software 
package.  For  additional  details  about  the  PRIME  software,  see  [21]. 

For  each  of  the  four  cases  of  simulated  data  detailed  in  Section  2.3,  two  interactive  Bezier 
distribution  fits  were  performed.  The  first  fit  for  each  case  uses  6  control  points  (requiring  the 
estimation  of  10  parameters)  while  the  second  fit  uses  12  control  points  (requiring  22  estimated 
parameters).  Figures  15  and  16  show  the  Bezier  fits  for  the  sodium/rectangular  case  using  6  control 
points  and  12  control  points,  respectively.  The  fits  in  these  figures  are  indicative  of  the  fits  obtained 
for  the  other  three  simulated  cases.  The  K-S  goodness-of-fit  statistic  is  0.01996  for  the  fit  with  6 
control  points  and  0.01535  for  the  fit  with  12  control  points.  Compared  to  the  K-S  goodness-of-fit 
values  for  the  corresponding  Johnson  fits  in  Tables  3  and  4,  the  Bezier  fits  appear  to  provide  a 
more  accurate  fit  to  the  simulated  data  since  the  maximum  discrepancy  between  the  empirical  and 
fitted  CDFs  is  much  smaller  for  the  Bezier  fits. 

From  Figure  16  it  is  evident  that  if  enough  control  points  are  used,  it  is  possible  to  fit  the 
simulated  (empirical)  data  very  closely  using  a  Bezier  function.  At  the  same  time,  however,  it  may 
not  be  desirable  to  track  the  empirical  data  exactly.  For  example,  the  empirical  data  in  Figure  16 
suggests  there  may  be  a  second  mode  in  the  lower  tail  of  the  distribution.  If  enough  control  points 
are  used,  it  is  possible  to  capture  the  bimodal  behavior  of  the  data  using  a  Bezier  density  function. 
Physically,  however,  there  is  not  any  evidence  to  suggest  that  the  underlying  true  density  of  cluster 
communication  distances  is  bimodal.  As  a  result,  it  is  not  clear  that  the  Bezier  fits  provide  a  better 
estimate  of  P(r)  than  the  Johnson  fits  in  Figures  7-14. 

The  primary  disadvantage  of  using  the  Bezier  method  is  the  property  that  the  corresponding 
PDF  (7)  has  a  much  more  complicated  form  than  either  the  cubic  spline  or  the  Johnson  PDF.  At 
this  point,  obtaining  stiffness  predictions  using  estimated  Bezier  density  functions  is  the  subject  of 
on-going  research. 
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(a)  (b) 


Figure  15:  Estimated  PDF  (a),  and  CDF  (b)  for  the  so dium/rect angular  case  using  a  Bezier 
function  with  6  control  points. 


(a)  (b) 


Figure  16:  Estimated  PDF  (a),  and  CDF  (b)  for  the  so  dium/rect  angular  case  using  a  Bezier 
function  with  12  control  points. 

4.  Concluding  Remarks 

In  this  paper  we  present  a  method  for  simulating  the  atomic  structure  of  the  Nafion  polymer  back¬ 
bone  chain  with  the  purpose  of  estimating  the  probability  density  function  P(r)  of  polymer  chain 
length.  Estimation  of  P(r)  is  the  first  (nanoscale)  step  in  a  multiscale  modeling  approach  for  the 
prediction  of  material  stiffness.  After  developing  an  appropriate  expression  for  P(r ),  the  resulting 
function  may  be  used  in  a  nracroscopic-level  model  to  compute  a  Young’s  modulus  through  applica¬ 
tion  of  Boltzmann’s  approach  to  statistical  thermodynamics  —  see  [22].  Successful  implementation 
of  P(r)  into  such  a  macroscopic  model  represents  a  significant  expansion  of  the  capabilities  of  RIS 
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theory.  In  order  to  achieve  this  goal  an  expression  for  P(r )  which  is  both  accurate  and  easily 
manipulated  mathematically  is  required. 

Thus,  in  addition  to  proposing  an  appropriate  simulation  for  multiscale  study  of  Nafion  stiffness, 
in  this  paper  we  have  also  proposed  two  methods  for  estimating  the  probability  density  function 
P(r)  as  alternatives  to  the  cubic  spline  method  used  in  [12].  Both  the  proposed  Johnson  and  Bezier 
density  estimation  methods  provide  a  more  flexible  approach  to  density  estimation  than  the  cubic 
spline  method.  Furthermore,  while  using  Bezier  curves  to  estimate  P(r)  offers  a  highly  flexible 
technique  for  estimating  a  wide  variety  of  distribution  shapes,  the  Johnson  family  of  distributions 
have  a  simple  four-parameter  functional  form.  Additionally,  the  parameters  of  the  Johnson  dis¬ 
tribution  have  a  specific  physical  meaning  and  current  research  efforts  are  focusing  on  identifying 
relationships  between  the  Johnson  parameters  and  predicted  stiffness. 
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